
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include "helper_timer.h"

#include <stdlib.h>
#include <stdio.h>

#include "util/mgpucontext.h"		// Include all MGPU kernels.
#include "moderngpu.cuh"		// Include all MGPU kernels.

#include "include/kernel.h"


#define COMPUTE_CLUSTERS

cudaError_t addWithCuda(int *c, const int *a, const int *b, unsigned int size);
cudaError_t ComputePersonalizedScoresRatesGPU(int userThreadsPerBlock);


namespace mgpu {
	////////////////////////////////////////////////////////////////////////////////
	// Device-side comparison operators

	struct lessFloat2 : public std::binary_function < float2, float2, bool > {
		MGPU_HOST_DEVICE bool operator()(float2 a, float2 b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
	};
}

namespace tests {

	/*
	v - Indicators(m x n) , m = companies, n = indicators no,  column ordered
	p  - weights (m x n) column ordered

	s - total score (1xm) , s[i]  = 0 if company i not in reference universe.
	*/

	void
		companiesScoresCPU(const float *v, const float *p, float * s, float* esg_score, int companiesNo, int indicatorsNo)
	{

		StopWatchInterface *timer = NULL;
		sdkCreateTimer(&timer);
		sdkResetTimer(&timer);

		sdkStartTimer(&timer);

		//47
		const int clusterStarts[52] = { 9, 10, 15, 16, 20, 29, 30, 32, 33, 34, 38, 43, 46, 48, 49, 55, 56, 57, 58, 59, 66, 67, 72, 73, 74, 87, 88, 89, 91, 92, 93, 94, 98, 99, 102, 107, 108, 110, 114, 117, 118, 125, 127, 129, 130, 135, 137, 138, 140, 141, 143, 147 }; // clusters interval (48 intervals)
		const int clusterGroups[53] = { 1, 0, 1, 2, 0, 1, 2, 1, 0, 2, 1, 3, 1, 3, 2, 1, 0, 1, 3, 2, 1, 3, 1, 2, 3, 1, 3, 2, 1, 3, 1, 0, 1, 0, 3, 1, 3, 0, 1, 3, 2, 1, 3, 1, 2, 3, 1, 3, 1, 3, 1, 3, 2 }; // group for each interval 0 - num groups <= intervals-1
		const int esgEnds[] = { 84, 100 };


		for (int line = 0; line < companiesNo; line++)
		{

			float score = 0;


			float clusterGroupsResults[4] = { 0, 0, 0, 0 };
			float clusterGroupScore = 0;
			int currentClusterGroup = 0;



			float groupRelativeWeight = 0;
			float groupScore = 0;
			char currentGroup = 0;

			int index = line;
			for (int i = 0; i < indicatorsNo; i++)
			{
				if (i == esgEnds[currentGroup])
				{
					esg_score[line + currentGroup * companiesNo] = groupScore / groupRelativeWeight;
					score += groupScore;
					++currentGroup;
					groupRelativeWeight = 0;
					groupScore = 0;
				}

				if (i == clusterStarts[currentClusterGroup])
				{
					int currentClusterIndex = clusterGroups[currentClusterGroup];
					clusterGroupsResults[currentClusterIndex] = clusterGroupScore + clusterGroupsResults[currentClusterIndex];
					++currentClusterGroup;
					clusterGroupScore = 0;

				}


				float pCurrent = p[index];
				float indicatorRating = pCurrent * v[index];

				groupRelativeWeight += pCurrent;
				groupScore += indicatorRating;

				clusterGroupScore += indicatorRating;
				//clusterGroupScore += 1;

				index = index + companiesNo;
			}

			esg_score[line + currentGroup * companiesNo] = groupScore / groupRelativeWeight;


			int currentClusterIndex = clusterGroups[currentClusterGroup];
			clusterGroupsResults[currentClusterIndex] = clusterGroupScore + clusterGroupsResults[currentClusterIndex];

			s[line] = score;
		}

		sdkStopTimer(&timer);


		printf("done.\n");
		printf("time spent executing by CPU : %.2f ms\n", sdkGetTimerValue(&timer));
	}

	__global__  void
		companiesScoresKernel(const float *v, const float *p, float * s, float* esg_score, int companiesNo, int indicatorsNo)
	{

		int line = (blockDim.x * blockIdx.x + threadIdx.x);
		if (line >= companiesNo)
			return;

		float score = 0;
		const int esgEnds[] = { 84, 100 };

#ifdef COMPUTE_CLUSTERS

		//47
		const int clusterStarts[52] = { 9, 10, 15, 16, 20, 29, 30, 32, 33, 34, 38, 43, 46, 48, 49, 55, 56, 57, 58, 59, 66, 67, 72, 73, 74, 87, 88, 89, 91, 92, 93, 94, 98, 99, 102, 107, 108, 110, 114, 117, 118, 125, 127, 129, 130, 135, 137, 138, 140, 141, 143, 147 }; // clusters interval (48 intervals)
		const int clusterGroups[53] = { 1, 0, 1, 2, 0, 1, 2, 1, 0, 2, 1, 3, 1, 3, 2, 1, 0, 1, 3, 2, 1, 3, 1, 2, 3, 1, 3, 2, 1, 3, 1, 0, 1, 0, 3, 1, 3, 0, 1, 3, 2, 1, 3, 1, 2, 3, 1, 3, 1, 3, 1, 3, 2 }; // group for each interval 0 - num groups <= intervals-1
		float clusterGroupsResults[4] = { 0, 0, 0, 0 }; // num group
		float clusterGroupScore = 0;
		int currentClusterGroup = 0;

#endif


		float groupRelativeWeight = 0;
		float groupScore = 0;
		char currentGroup = 0;

		int index = line;
		for (int i = 0; i < indicatorsNo; i++)
		{
			if (i == esgEnds[currentGroup])
			{
				esg_score[line + currentGroup * companiesNo] = groupScore / groupRelativeWeight;
				score += groupScore;
				++currentGroup;
				groupRelativeWeight = 0;
				groupScore = 0;
			}

#ifdef COMPUTE_CLUSTERS
			if (i == clusterStarts[currentClusterGroup])
			{
				//clusterGroupsResults[currentClusterIndex] = clusterGroupScore;
				clusterGroupsResults[clusterGroups[currentClusterGroup]] = clusterGroupScore + clusterGroupsResults[clusterGroups[currentClusterGroup]];
				++currentClusterGroup;
				clusterGroupScore = 0;

			}
#endif

			float pCurrent = p[index];
			float indicatorRating = pCurrent * v[index];

			groupRelativeWeight += pCurrent;
			groupScore += indicatorRating;

#ifdef COMPUTE_CLUSTERS
			clusterGroupScore += indicatorRating;
			//clusterGroupScore += 1;
#endif

			index = index + companiesNo;
		}


		esg_score[line + currentGroup * companiesNo] = groupScore / groupRelativeWeight;

#ifdef COMPUTE_CLUSTERS
		int currentClusterIndex = clusterGroups[currentClusterGroup];
		clusterGroupsResults[currentClusterIndex] = clusterGroupScore;
#endif

		//s[line] = score;

		s[line] = companiesNo - line;
	}

	__global__  void
		companiesScoresKernelEsgInline(const float *v, const float *p, float * s, float* esg_score, int companiesNo, int indicatorsNo)
	{

		int line = (blockDim.x * blockIdx.x + threadIdx.x);
		if (line >= companiesNo)
			return;

		float score = 0;

		const int starts[] = { 84, 100 };
		float clusterResults[3] = { 0, 0, 0 }; // num group
		//const int g_start = 100;

		float groupRelativeWeight = 0;
		float groupScore = 0;


		int index = line;
		for (int i = 0; i < starts[0]; i++)
		{
			groupScore += p[index] * v[index];
			groupRelativeWeight += p[index];
			index = index + companiesNo;
		}

		clusterResults[0] = groupScore / groupRelativeWeight;
		score += groupScore;
		groupRelativeWeight = 0;
		groupScore = 0;


		for (int i = starts[0]; i < starts[1]; i++)
		{
			groupScore += p[index] * v[index];
			groupRelativeWeight += p[index];
			index = index + companiesNo;
		}

		clusterResults[1] = groupScore / groupRelativeWeight;

		score += groupScore;
		groupRelativeWeight = 0;
		groupScore = 0;

		for (int i = starts[1]; i < indicatorsNo; i++)
		{
			groupScore += p[index] * v[index];
			groupRelativeWeight += p[index];
			index = index + companiesNo;
		}

		clusterResults[2] = groupScore / groupRelativeWeight;

		esg_score[line] = clusterResults[0];
		esg_score[line + companiesNo] = clusterResults[1];
		esg_score[line + 2 * companiesNo] = groupScore / groupRelativeWeight;

		score += groupScore;
		s[line] = score;
	}

	__global__  void
		companiesScoresKernelClusters(const float *v, const float *p, int companiesNo, int indicatorsNo)
	{

		int line = (blockDim.x * blockIdx.x + threadIdx.x);
		if (line >= companiesNo)
			return;

		float score = 0;

		const int clusterStarts[] = { 12, 68 }; // clusters interval  0-12-68-end (3 intervals)
		const int clusterGroups[3] = { 0, 1, 0 }; // group for each interval 0 - num groups <= intervals-1
		float clusterGroupsResults[2] = { 0, 0 }; // num group


		float clusterGroupScore = 0;
		int currentClusterGroup = 0;

		for (int i = 0; i < indicatorsNo; i++)
		{

			if (i == clusterStarts[currentClusterGroup])
			{
				int currentClusterIndex = clusterGroups[currentClusterGroup];
				clusterGroupsResults[currentClusterIndex] = clusterGroupScore;
				++currentClusterGroup;
				clusterGroupScore = clusterGroupsResults[currentClusterIndex];
			}

			int index = companiesNo * i + line;
			clusterGroupScore += p[index] * v[index];

		}


		int currentClusterIndex = clusterGroups[currentClusterGroup];
		clusterGroupsResults[currentClusterIndex] = clusterGroupScore;
	}

	/*
	u - reference universe U(1 x m), row ordered, 1 company is the univers, 0 company is not in the universe,

	filter : 2bits pt  + 2bits tt  private| public + type A|type B
	attributes: 4 bits: 2 pt, 2tt private| public + type A|type B
	pt, tt bits - flags.

	*/
	__global__  void
		filterCompaniesKernel(int companiesNo, const char * u, char* attributes, char filter, char *x)
	{
		int line = (blockDim.x * blockIdx.x + threadIdx.x);
		if (line >= companiesNo)
			return;

		x[line] = u[line] & attributes[line] & filter;
	}

	//s* - score per company, -1 daca nu indeplineste conditia. ordonat dupa grup.
	//N number of companies
	void SegmentedSort_Pairs(mgpu::CudaContext& context, float *s, int *companiesMappings, int *segmentsHeads, int N, int numSegs)
	{
		mgpu::SegSortPairsFromIndices(s, companiesMappings, N, segmentsHeads, numSegs, context);
	}

	cudaError_t ComputePersonalizedScoresRatesGPU(int userThreadsPerBlock)
	{
		mgpu::ContextPtr context = mgpu::CreateCudaDevice(0);

		cudaEvent_t start, stop;
		cudaError_t cudaStatus = cudaEventCreate(&start);
		if (cudaStatus != cudaSuccess)
			return cudaStatus;

		cudaStatus = cudaEventCreate(&stop);

		int companiesNo = 5128 * 1;
		int indicatorsNo = 148;
		int numElements = companiesNo * indicatorsNo;

		size_t size = numElements * sizeof(float);
		float *v_h = (float *)malloc(size); // indicators score
		float *p_h = (float *)malloc(size); // relative Weights matrix

		float *s_h = (float *)malloc(companiesNo * sizeof(float)); // companies scores vector
		float *group_score_h = (float *)malloc(3 * companiesNo * sizeof(float));
		char *attributes_h = (char*)malloc(companiesNo * sizeof(char)); // companies scores vector
		char *u_h = (char*)malloc(companiesNo * sizeof(char)); // companies scores vector
		char *x_h = (char*)malloc(companiesNo * sizeof(char)); // companies scores vector

		float* esgScores[3];
		for (int i = 0; i < 3; i++)
		{
			float * group_score_h = (float *)malloc(companiesNo * sizeof(float)); // companies scores vector
		}


		printf("executing companies = %d  , indicators = %d ... .\n", companiesNo, indicatorsNo);

		for (int i = 0; i < companiesNo; ++i)
		{
			attributes_h[i] = 1; // 1 | 2 + 4| 8
			u_h[i] = i % 2;
			for (int j = 0; j < indicatorsNo; ++j)
			{
				int currentIndex = indicatorsNo * i + j;
				p_h[currentIndex] = rand() / (float)RAND_MAX + 1;
				v_h[currentIndex] = rand() / (float)RAND_MAX + 1;
			}
		}

		float *v_d = NULL; //indicators matrix
		cudaStatus = cudaMalloc((void **)&v_d, size);
		if (cudaStatus != cudaSuccess)
			return cudaStatus;

		cudaStatus = cudaMemcpy(v_d, v_h, size, cudaMemcpyHostToDevice);
		if (cudaStatus != cudaSuccess)
			return cudaStatus;

		float *p_d = NULL; //relative weights matrix
		cudaStatus = cudaMalloc((void **)&p_d, size);
		if (cudaStatus != cudaSuccess)
			return cudaStatus;

		cudaStatus = cudaMemcpy(p_d, p_h, size, cudaMemcpyHostToDevice);
		if (cudaStatus != cudaSuccess)
			return cudaStatus;

		char *attributes_d = NULL; //relative weights matrix
		cudaStatus = cudaMalloc((void **)&attributes_d, companiesNo* sizeof(char));
		cudaStatus = cudaMemcpy(attributes_d, attributes_h, companiesNo * sizeof(char), cudaMemcpyHostToDevice);

		char *u_d = NULL; //relative weights matrix
		cudaStatus = cudaMalloc((void **)&u_d, companiesNo* sizeof(char));
		cudaStatus = cudaMemcpy(u_d, u_h, companiesNo * sizeof(char), cudaMemcpyHostToDevice);

		float *s_d = NULL;// companies scores vector
		cudaStatus = cudaMalloc((void **)&s_d, companiesNo * sizeof(float));

		char *x_d = NULL;// companies scores vector
		cudaStatus = cudaMalloc((void **)&x_d, companiesNo * sizeof(char));


		float *group_score_d = NULL;
		cudaStatus = cudaMalloc((void **)&group_score_d, 3 * companiesNo * sizeof(float));


		int threadsPerBlock = 32 * 6;
		if (userThreadsPerBlock > 0)
		{
			threadsPerBlock = userThreadsPerBlock;
		}

		int blocksPerGrid = companiesNo / threadsPerBlock;
		if (companiesNo % threadsPerBlock != 0)
		{
			++blocksPerGrid; // ensure that all companies assigned to a thread.
		}

		//start compute scores async

		StopWatchInterface *timer = NULL;
		sdkCreateTimer(&timer);
		sdkResetTimer(&timer);

		sdkStartTimer(&timer);

		///////////////-------------------  compute score ---------------------

		//filterCompaniesKernel(int companiesNo, const char * u, char* attributes, char filter, char *x)
		filterCompaniesKernel << <blocksPerGrid, threadsPerBlock >> >(companiesNo, u_d, attributes_d, 1, x_d);

		///companiesScoresKernel(const float *v, const float *p, float * s, float* esgScores[3], int companiesNo, int indicatorsNo)
		companiesScoresKernel << <blocksPerGrid, threadsPerBlock >> >(v_d, p_d, s_d, group_score_d, companiesNo, indicatorsNo);


		MGPU_MEM(int) s_companies_map = context->FillAscending<int>(companiesNo, 0, 1);
		const int groupsCount = 10;//companiesNo / 1000;
		const int SegHeads[groupsCount] = { 3, 100, 200, 500, 1000, 2000, 2100, 2200, 2300, 5000 };
		MGPU_MEM(int) companies_groups = context->Malloc(SegHeads, groupsCount);
		//MGPU_MEM(int) companies_groups = context->FillAscending<int>(groupsCount, 0, groupsCount);


		// Record the start event
		cudaStatus = cudaEventRecord(start, NULL);

		//// --- s sorted by score s_companies_mem
		SegmentedSort_Pairs(*context, s_d, s_companies_map->get(), companies_groups->get(), companiesNo, groupsCount);

		cudaStatus = cudaEventRecord(stop, NULL);
		cudaStatus = cudaEventSynchronize(stop);

		float msecTotal = 0.0f;
		cudaStatus = cudaEventElapsedTime(&msecTotal, start, stop);

		cudaStatus = cudaMemcpy(s_h, s_d, companiesNo * sizeof(float), cudaMemcpyDeviceToHost);
		sdkStopTimer(&timer);
		printf("done.\n");
		printf("CPU : %.2f ms\n", sdkGetTimerValue(&timer));
		printf("Time= %.3f msec  threadsPerBlock= %d  BlocksPerGrid=%d\n", msecTotal, threadsPerBlock, blocksPerGrid);


		printf("Executiong on cpu: ... \n");
		companiesScoresCPU(v_h, p_h, s_h, group_score_h, companiesNo, indicatorsNo);


		free(v_h);
		free(p_h);
		free(s_h);

		cudaFree(v_d);
		cudaFree(p_d);
		cudaFree(s_d);
	}


	template<typename T>
	void BenchmarkSortKeys(int count, int numIt, mgpu::CudaContext& context) {
#ifdef _DEBUG
		numIt = 1;
#endif

		MGPU_MEM(T) source = context.GenRandom<T>(count, 0, (T)count);
		MGPU_MEM(T) data = context.Malloc<T>(count);
		std::vector<T> sourceHost;
		source->ToHost(sourceHost);

		double mgpuElapsed = 0;
		for (int it = 0; it < numIt; ++it) {
			source->ToDevice(data->get(), count);
			context.Start();
			MergesortKeys(data->get(), count, mgpu::less<T>(), context);
			mgpuElapsed += context.Split();
		}

		std::stable_sort(sourceHost.begin(), sourceHost.end());
		double cpuElapsed = context.Split();

		double bytes = sizeof(T) * count;
		double mgpuThroughput = (double)count * numIt / mgpuElapsed;
		double mgpuBandwidth = bytes * numIt / mgpuElapsed;
		double cpuThroughput = count / cpuElapsed;
		double cpuBandwidth = bytes / cpuElapsed;

		printf("%s: %9.3lf M/s  %7.3lf GB/s   %9.3lf M/s  %7.3lf GB/s\n",
			mgpu::FormatInteger(count).c_str(), mgpuThroughput / 1.0e6,
			mgpuBandwidth / 1.0e9, cpuThroughput / 1.0e6,
			cpuBandwidth / 1.0e9);

		printf("\t Elapsed ms  Gpu:%f    CPU:%f \n\n",
			mgpuElapsed * 1000,
			cpuElapsed * 1000
			);


		// Verify
		std::vector<T> host;
		data->ToHost(host);
		for (int i = 0; i < count; ++i)
			if (sourceHost[i] != host[i]) {
				printf("MISMATCH AT %d\n", i);
					exit(0);
			}
	}

	void
		ComputeMgpuSort(int argc, char** argv){
		mgpu::ContextPtr context = mgpu::CreateCudaDevice(argc, argv, true);

		int N = 512;

		for (int i = 1; i <= 5; i++){

			MGPU_MEM(float2) vals = context->Malloc<float2>(N);
			std::vector<float2> sourceHost;
			vals->ToHost(sourceHost);

			context->Start();
			mgpu::MergesortKeys<float2>(vals->get(), N, mgpu::lessFloat2(), *context);
			double timeS = context->Split();


			std::stable_sort(sourceHost.begin(), sourceHost.end(), mgpu::lessFloat2());
			double cpuElapsed = context->Split();

			printf("\nSorted values in : %fms  & cpu: %fms\n", timeS * 1000, cpuElapsed * 1000);

			BenchmarkSortKeys<float>(N, 1, *context);
			N = N * 10;
			//vals->Release();
		}
	}
}

main(int argc, char** argv)
{
//	tests:ComputePersonalizedScoresRatesGPU(0);


	cudaEsg::ComputeContext computeContext;
	computeContext.InitializeWithRandomData();
	computeContext.Execute(0, 0, 127); //no filter first univers

	cudaDeviceReset();
	//ComputeMgpuSort(argc, argv);

	//int N = 10;
}

